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We present long-term-stable and convergent evolntions of head-on black hole collisions and ex¬ 
traction of gravitational waves generated dnring the merger and subseqnent ring-down. The new 
ingredients in this work are the use of fixed mesh-refinement and dynamical singularity excision tech¬ 
niques. We are able to carry out head-on collisions with large initial separations and demonstrate 
that our excision infrastructure is capable of accommodating the motion of the individual black 
holes across the computational domain as well as their merger. We extract gravitational waves from 
these simulations using the Zerilli-Moncrief formalism and find the ring-down radiation to be, as 
expected, dominated by the £ = 2, m = 0 quasi-normal mode. The total radiated energy is about 
0.1 % of the total ADM mass of the system. 

PACS numbers: 04.25Dm, 04.30Db 

I. INTRODUCTION 

One of the most urgent problems currently under in¬ 
vestigation in general relativity is the coalescence of bi¬ 
nary compact objects. The importance of this scenario 
largely arises from its significance for the new research 
field of gravitational wave astronomy; ground-based grav¬ 
itational wave interferometric detectors, such as LIGO, 

GEO600, VIRGO and TAMA, have started collecting 
data. The in-spiral and merger of neutron stars and black 
holes is considered among the most promising sources of 
gravitational radiation to be detected by these instru¬ 
ments. Theoretical predictions of the resulting wave pat¬ 
terns will be important both for the detection and for 
the physical interpretation of such measurements. Al¬ 
though approximation techniques have been remarkably 
successful, the theoretical calculation of wave forms from 
compact object binaries will most likely require the nu¬ 
merical solution of the Einstein held equations. 

Very important milestones have been achieved recently 
in the numerical modeling of compact object binaries. 

For instance, a few orbits of binary neutron stars are 
now possible (see e.g. HUlli) Also, with binary black 
holes, Briigmann et al. Q have been able to carry out 
evolutions lasting for timescales longer than one orbital 
period. These varied efforts in modeling compact binary 
systems with numerical relativity share a common goal 
- to extract accurate waveform templates (see @ for a 
study on accuracy requirements of such waveform calcu¬ 
lations). This “Holy Grail” of numerical relativity, how¬ 
ever, still eludes the community. A unique problem faced 
in the simulation of black hole spacetimes is the presence 
of physical singularities, which makes this task particu¬ 


larly delicate from a numerical point of view. 

For a long time after the pioneering work of Eppley 
and Smarr in the 1970s idi , progress in numeri¬ 
cal black hole evolutions was rather slow, and the codes 
suffered from instabilities after relatively short times. 
In contrast, the last few years have seen a change of 
strategy in the area of numerical relativity. More em¬ 
phasis has been put on analyzing the structure of the 
evolution equations und erlYing the numerical codes (see 
e.g. in combination 

with increased computational resources and advanced 
treatments of the spacetime singularities, this has led to 
vastly improved stability properties in the evolution of 
spacetimes with a single black hole either at rest or mov¬ 
ing through the computational domain m m m 113 . 
In view of these encouraging results, attention has been 
shifting again toward the application of these techniques 
to the evolution and merger of binary black holes (see 

eg. mu mill El). 

In this paper, we will focus on the numerical evolution 
of head-on collisions of black holes implemented in the 
framework of dynamical singularity excision. The study 
of head-on colliding black holes has attracted a good deal 
of attention in the past, both numerically and analyt¬ 
ically. We have already mentioned the early work of 
Eppley and Smarr who numerically evolved Misner data 
for equal-mass black hole collisions in axisymmetry. The 
waveforms extracted from these evolutions differ little in 
amplitude and phase from those obtained in perturbative 
calculations |3 . The head-on collision was re-investigated 
numerically in the 1990s El El El El E3 with gener¬ 
ally good agreement between numerical results and those 
obtained from perturbation theory for various initial con- 
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figurations. As a whole, the approximate analytic study 
of binary black hole head-on collisions has been remark¬ 
ably successful and has generated a wealth of literature. 
The late stage of the collision, when the separation of 
the two holes is sufficiently small, is well described by 
the close limit approximation (see e.g. 1 ^ as well 
as [Sli) for an estimate of the quality of the approxima¬ 
tion using second-order perturbation theory). Similarly 
the particle approximation has been used to describe the 
in-fall phase with large black-hole separation either in 
the framework of post-Newtonian theory (see e.g. |^i 
or using Newtonian trajectories with appropriate modi¬ 
fications In particular, Anninos et al. com¬ 

bined their “particle-membrane” technique with results 
from close limit calculations to provide estimates of the 
total radiated energy that are in good agreement with 
the results of numerical relativity over the whole range 
of initial black hole separations (see their Fig. 11). 

The encouraging results of these perturbative approxi¬ 
mations gave rise to the idea of combining them with full 
numerical evolutions via matching on slices of constant 
time. With this temporal matching, it was possible to 
extend the life of simulations of binary black holes which 
would otherwise stop because of numerical instabilities. 
This is now commonly known as the “Lazarus” approach; 
developed by Baker et al. |3^ it has demonstrated good 
waveform agreement with non-linear evolutions. 

The general picture that has emerged in analytic as 
well as numerical studies of equal mass head-on collisions 
is that gravitational wave emission is dominated by the 
(. = 2^ m = Q quasi-normal oscillations of the post-merger 
black hole and that the total radiated wave energy is 
typically less than a percent of the system’s ADM mass. 

In spite of their limited astrophysical relevance, simu¬ 
lations of black hole head-on collisions still represent an 
important problem. In view of the eventual target, the 
in-spiral and merger of a binary black hole, these sim¬ 
ulations constitute a key milestone in the development 
of numerical codes, a natural next step after successful 
evolutions of single black holes. While numerical codes 
have continued to improve in the simulation of head-on 
or grazing collisions |3a 113 , long-term-stable simu¬ 

lations have only been obtained rather recently by Al- 
cubierre et al. |23 using the puncture method. In 
the combination of this approach with the “Simple Ex¬ 
cision” technique of Alcubierre and Briigmann |19| has 
been shown to preserve the long-term stability and yield 
good agreement in the wave forms. Almost parallel to 
this work black hole head-on collisions have been sim¬ 
ulated in the framework of mesh refinement and 
fourth-order finite differencing techniques |39|. In this 
paper we present the first such long-term-stable evolu¬ 
tions obtained in fixed mesh-refinement as well as using 
the dynamical singularity excision technique. 

The simulations presented in this paper have been ob¬ 
tained with the Maya evolution code. This code and 
in particular its dynamical sing ularity excision technique 
are described in detail in (referred to as Paper 1 


from now on). A key development in recent months has 
been the combination of the Maya code with the mesh- 
refinement package Carpet ^3- Carpet has enabled 
us to push the outer boundaries substantially further out 
than was possible with the uni-grid version of Maya; this 
was crucial in the extraction of gravitational wave signals. 

The paper is organized as follows: The construction 
of binary black hole initial data used in this work is 
described in Sec. im while the details of the numerical 
evolution scheme are given in Sec. cni The gravita¬ 
tional waveforms are analyzed in Sec. nYi with regard to 
their convergence properties and the physical informa¬ 
tion they provide. We conclude in Sec. m with a discus¬ 
sion of our results and the implications for future work. 
In the appendix we present the details of the wave ex¬ 
traction using the Zerilli-Moncrief formalism on a sin¬ 
gle black hole background in Kerr-Schild coordinates. 
Throughout the paper, units are given in terms of M, 
with M = + ^m, the sum of the bare masses of the 

individual black holes in the binary. 

II. INITIAL DATA 

The construction of astrophysically meaningful initial 
data for simulations of compact binaries is a highly non¬ 
trivial problem and represents an entire branch of re¬ 
search in numerical relativity (see for a detailed 
overview). 

For the simulations discussed in this work, we con¬ 
struct initial data according to the procedure of Matzner 
et al. ^3 j referred to from now on as “HuMaSh” data. 
This type of data has been used, for example, in binary 
grazing-collision evolutions by Brandt et al. |^. The 
starting point of the HuMaSh data is the Kerr-Schild 
spacetime metric 


9tJ.v 

— “ 1 “ 2 ^ ; 

( 1 ) 

where rj^v is the Minkowski metric, and L 

is a null-vector 

with respect to both and 77 ^ 1 ,. In terms of 3-1-1 quan- 

titles, the Kerr-Schild metric is given by 
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In particular, for a non-rotating black hole of mass m at 

rest at the origin 
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with r = Sij xi. 
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As pointed out in Ref. Kerr-Schild coordinates 
have two attractive features. First, with relevance to 
singularity excision, the coordinates penetrate the hori¬ 
zon. Second, the Kerr-Schild spacetime metric is form- 
invariant under a boost transformation. Specifically, un¬ 
der a boost V along the z-axis, all that changes is 


^{R — vZ) X y j {Z — V R) 

^ ^ [ R ’R’R' R J’ 

where 7 = 1/a/1 — Z = 7(2 — vt), and R"^ = + 

J/2 + 

A spatial metric approximating initial data of two 
black holes, with bare masses and respectively, 

can be given by 

lij ~ lij ~ ^ij ) ( 8 ) 

where the "^ 7 ^ and '^ 7 ^- are to be understood as the 
Kerr-Schild spatial metric @ evaluated at the corre¬ 
sponding black hole coordinate position and boost ve¬ 
locity. Following the HuMaSh prescription, the extrinsic 
curvature is constructed from 

, (9) 

where the j and ^ j are also to be understood 
as the Kerr-Schild extrinsic curvature evaluated at the 
corresponding black hole coordinate position and boost 
velocity. 

It is important to point out that the HuMaSh data 
do not provide a close limit; that is, as the coordinate 
separation of the black holes shrinks to zero, the resulting 
metric is not that of a single black hole with mass M = 
+ ^m. Furthermore, as given by Q and @, the 
HuMaSh data do not satisfy the Einstein constraints. As 
suggested by Marronetti et al. could attenuate 

the constraint violations of the HuMaSh data using an 
exponential damping factor. Bonning et al. used 

this attenuating procedure and provided a fairly detailed 
investigation of the general properties of the underlying 
HuMaSh data. In particular, they show that it has the 
correct gravitational binding energy in the Newtonian 
far-separation limit. 

By construction the HuMaSh data will satisfy the con¬ 
straint equations only in the limit of infinite separation 
of the two black holes. A rigorous treatment of a sce¬ 
nario with finite separation would thus require using the 
HuMaSh data as freely-specifiable background data and 
then solving the Hamiltonian and momentum constraint 
according to York’s initial data prescription Since 
one of our primary goals is to extract gravitational waves 
from the black hole simulations, we need a sufficiently 
large computational domain, only attainable via mesh- 
refinement. Unfortunately, elliptic solvers operating on 
numerical grids with various refinement levels have not 
yet been developed in our Maya code. Using the uni-grid 
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FIG. 1: Convergence analysis of the Hamiltonian constraint 
along the z-axis at t = 0 (upper panel) and t = 25.92 M 
(lower panel). 


solvers available, on the other hand, is computationally 
prohibitive. Therefore, as with the grazing collision work 
in Ref. |23| . the simulations in this paper were obtained 
using the constraint-violating HuMaSh initial data. 

In this context one has to keep in mind that the con¬ 
straints will always be violated at the level of numerical 
accuracy in a real simulation. In consequence, the con¬ 
straint violations inherent to the HuMaSh construction 
are problematic only if they are significant relative to 
the limit set by finite grid resolutions in the course of 
the evolution. Furthermore, the constraint violations are 
reduced as we increase the initial separation of the black 
holes. 

In order to estimate the significance of the constraint 
violations in the HuMaSh data relative to the numeri¬ 
cal evolution errors, we have analyzed the convergence 
properties of the constraint violations in the course of 
the simulations. Because the initial data are not exactly 
constraint satisfying, we do not expect the values of the 
constraints to converge to zero at second order in the con¬ 
tinuum limit. Let TLq denote the Hamiltonian constraint 
violation at the continuum level inherent to the HuMaSh 
data. Tfo is thus independent of the grid resolution. For 
a second-order convergent evolution code, the numerical 
constraint violation T-Lh will be given by 

nh=Ho + ch^, ( 10 ) 

where h denotes the grid resolution and c is a function 
of space and time, independent of h. We have carried 
out two head-on collisions with black holes separated by 
a coordinate distance 6 M. These runs have resolutions 
hi = 0.135M and /i 2 = hi/l.'bM on the finest refine¬ 
ment level respectively, with the grid spacing doubled in 
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about z = 22 M. We emphasize that these outer bound¬ 
ary effects are not related to the question of initial data 
and have also been observed, for example, in the punc¬ 
ture evolutions of Alcubierre et al. . This viewpoint is 
further strengthened by the convergence analysis of the 
waveforms in Sec. lTV Rl which also reveal the propagation 
of noise from the outer boundaries as the limiting factor. 

In order to shed more light on the constraint violations 
at small z, we show in Fig.[21the Hamiltonian constraint 
violations along the a;-axis at the same times. Due to 
the axisymmetry in the problem, the results along the 
y-axis are identical. Clearly the initial data do not ex¬ 
hibit second-order convergence, in agreement with the 
results of Fig.^ near z = 0. In the course of the evo¬ 
lution, however, as the discretization error increases and 
dominates the total violations we do find second-order 
convergence of the Hamiltonian constraint except for the 
outer boundary effects mentioned above. 


FIG. 2: Convergence analysis of the Hamiltonian constraint 
along the ®-axis at t = 0 (upper panel) and t = 25.92 M 
(lower panel). 


each of the three additional coarser levels. (These are the 
“coarse” and “medium” simulations in the convergence 
analysis of the waveforms in Sec. ME) For a second- 
order convergent code we expect the constraint violations 
at the two resolutions to obey the relation 


Hhi _ 1-10+ ch\ 
^112 Ho + chi 


( 11 ) 


It is instructive to consider the following two limiting 
cases of this equation. First, if the contribution Ho due 
to the use of initially unsolved data dominates the total 
constraint violation, we obtain q — 1 and should observe 
no convergence at all. On the other hand, if Ho is negli¬ 
gible, we expect second-order convergence and q = 1.5^. 


In view of this result, we have investigated the point- 
wise convergence of both the Hamiltonian and the mo¬ 
mentum constraints. Fig.^ shows the Hamiltonian con¬ 
straint violation Hh along the z-axis (collision axis) for 
both runs, with the medium resolution result (/12 = 
hi/1.5) scaled by a factor of 1.5^. The upper panel 
shows Hh at t = 0 and the lower panel at t = 25.92 M. 
It is clear from this figure that at t = 0 second-order 
convergence is obtained in the strong field area close to 
the hole as well as in the far-field limit. At intermedi¬ 
ate distances, on the other hand, the convergence dete¬ 
riorates, in particular near z = 0, the region between 
the two holes. In the course of the evolution, however, 
the convergence approaches second order everywhere ex¬ 
cept in the domain of influence of the outer boundary, 
as clearly demonstrated in the lower panel of Fig.Q] The 
constraint violations are second-order convergent out to 


The results obtained for the momentum constraints 
show the same behavior, though they exhibit a stronger 
level of noise emanating from the refinement boundaries. 
We currently do not have an explanation for this discrep¬ 
ancy between Hamiltonian and momentum constraint 
but emphasize that the overall convergence of the mo¬ 
mentum constraint is still close to second order up to the 
outer boundary effects. 

Similar convergence properties are found for the £ 2 - 
norms of the Hamiltonian constraint H, the momen¬ 
tum constraints Mi and the auxiliary constraint = 
r* — In Fig.^lwe show the resulting plots for the 

x-components Mx and . As before, the constraint vio¬ 
lations obtained with fine resolution /12 are scaled by the 
expected factor 1.5^. The results for the other compo¬ 
nents of the constraint variables are similar and overall 
we find convergence close to second order after a tran¬ 
sition time of about 10 M. Outer boundary effects do 
not manifest themselves as clearly in the £ 2 -norms as in 
the point-wise convergence analysis. We attribute this 
feature to the fact that the norms are dominated by the 
large constraint violations near the black holes. 

In summary, our results confirm that the total con¬ 
straint violation is indeed dominated by the numerical 
inaccuracies inevitably generated during time evolutions 
at currently available grid resolutions. This analysis gen¬ 
eralizes the studies of Marronetti et al. to include 
the discretization errors arising from the time evolution; 
these lead to an extension of the allowed range of grid 
spacings from their limiting value of M/4 to about M/8 
as used in this work. While constraint solving of the 
initial data will be desirable for future simulations with 
higher resolutions, its effect on the overall constraint vio¬ 
lations encountered in our simulations would be marginal 
at most. We emphasize in this context that this is only 
the case because of the absence of exponentially growing 
constraint violating modes in our simulations. 
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FIG. 3: Convergence analysis of the ^ 2 -nornis as functions 
of time of the a;-component of the momentum and auxiliary 
constraints. 


III. NUMERICAL EVOLUTION 

From the operational point of view, it is convenient 
to divide a black hole head-on collision into pre-merger, 
merger and post-merger stages. In the pre-merger stage, 
the evolution contains two separate apparent horizons 
and therefore separate excision regions. The merger stage 
is signaled by the appearance of a single, highly dis¬ 
torted apparent horizon. Initially, these deformations 
of the horizon prevent the merger of the excision masks 
into a single spherical excision. The end of the merger 
stage is marked by the merger of excision masks when 
the deformations of the common apparent horizon have 
damped-out. The main constraint on the construction 
of the merged mask is that it must be completely inside 
the apparent horizon and fully contain the previous two 
disjoint excision regions. During the final stage (post¬ 
merger), the evolution consists of the settling or ring- 
down of the resulting black hole. 

There are certain aspects besides the handling of the 
excision region, such as gauge conditions, that require a 
different treatment during each of these stages. Before we 
discuss these differences in detail, we address those ingre¬ 
dients of the numerical evolution that remain unchanged 
throughout the complete simulation. 

The specific implementation of the 3-1-1 BSSN formu¬ 
lation and the black hole excision technique used in the 
Maya code have been described in detail in Paper 1. 
Here we provide a brief summary with particular regard 
to the changes necessary for working with fixed-mesh re¬ 
finement. 

We evolve the Einstein equations on a set of nested 
non-moving, uniform Cartesian grids with the resolution 


increasing by a factor of two for each additional refine¬ 
ment level. The Einstein field equations are evolved ac¬ 
cording to the BSSN 3-1-1 formulation |^0|. The ex¬ 
act shape of the equations as implemented in the code 
is given by Eqs. (6)-(10) in Paper 1; in particular this 
includes the additional term in the evolution of P® intro¬ 
duced by Yo et al. (the free parameter y present in 
that term is set to 2/3 for all simulations discussed in this 
paper). These equations are evolved in time with the it¬ 
erated Crank-Nicholson scheme. Spatial discretization is 
second-order centered finite differencing with the excep¬ 
tion of advection derivatives contained in the Lie deriva¬ 
tives; these are approximated with one-sided second- 
order stencils. With regard to the treatment of refine¬ 
ment boundaries, we follow the recipe of Schnetter et 
al. EH to avoid loss of convergence at these boundaries. 
Eor our code, this implies the use of a total of six buffer 
zones; two points (because of the one-sided derivatives 
mentioned above) for each of the three Crank-Nicholson 
iterations. 

With the exception of gauge quantities specified alge¬ 
braically as functions of the spacetime coordinates, we 
apply a Sommerfeld-like outgoing condition at the outer 
boundary . At the excision boundary, we extrapolate 
all evolution fields from the computational interior using 
quadratic or, when populating previously excised points, 
cubic polynomials. With regard to mesh refinement we 
note that the black hole excision is actively implemented 
on the finest level, only migrating to the coarser levels via 
the restriction of grid functions; consequently, we always 
require that the black holes remain within the range of 
the finest grid. In practice this is not a limitation since 
it is precisely the black hole neighborhood where steep 
gradients are present and high resolution is therefore re¬ 
quired. 

In Paper 1, the uni-grid version of this code has been 
shown to produce long-term-stable evolutions of single 
Schwarzschild black holes as well as evolutions of moving 
black holes lasting for about 130 M. In (referred to 
as Paper 2 from now on), we managed to extend the life 
times of these runs to at least 6000 M by describing the 
slicing condition in terms of a densitized lapse function 
[cf. Eqs. (6)-(8) therein] and demonstrated second-order 
convergence of the code. 

The remaining aspects of the numerical simulation are 
specific to the individual stages of the evolution and will 
now be discussed in turn. The main differences encoun¬ 
tered here have to do with gauge and excision. 


A. Pre-Merger Evolution 

Pre-merger, our spacetime contains two apparent hori¬ 
zons. Within each apparent horizon, there is a curvature 
singularity, hidden from us by singularity excision. As 
the holes approach each other, the excised regions - or 
excision masks - are allowed to move also. The new 
locations of the excision masks must be found by some- 
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how tracking the locations of the black hole singulari¬ 
ties. One possibility is to track the singularities via an 
apparent horizon finder; we have found, however, that 
a simple Gaussian fit to the shape of the trace of the 
extrinsic curvature is equally effective for this purpose. 
As the excision regions move, grid points that were pre¬ 
viously excised become part of the computational do¬ 
main. These uncovered points must be supplied with 
valid data for all fields; this “re-population” of points 
is achieved through extrapolation onto the mask surface 
along coordinate normals to the surface. Both the Gaus¬ 
sian mask tracking and the re-population procedure have 
been treated in detail in Papers 1 and 2. 

Besides the above, the current problem of evolving bi¬ 
nary data raises new challenges. Primary among these 
is the question of gauge conditions. While substan¬ 
tial progress has been made in the recent past in de¬ 
riving successful gauge conditions for evolutions that 
anchor black holes to a fixed coordinate location (see 
e- g- 0, E, 113; E3) EM EM) question as to suit¬ 
able gauge conditions for strongly time varying scenar¬ 
ios, such as spacetimes with moving black holes, remains 
largely unanswered. 

In view of the successful simulations of single mov¬ 
ing black holes in Paper 2, we have decided to pursue a 
similar approach for the binary system at hand, namely 
analytic shift and densitized lapse. A key benefit of using 
a densitized version of the lapse is that it facilitates an 
algebraic specification of the slicing condition, in terms of 
an analytic function of the spacetime coordinates, while 
satisfying criteria for strong and symmetric hyperbolic- 
ity as derived in or [^. Using a densitized lapse 
thus leads to a well-posed Gauchy problem as defined in 
[^. In contrast to the single black hole simulations of 
Paper 2, where algebraic forms of the shift and densi¬ 
tized lapse were available from the analytic single black 
hole solution, gauge conditions for our head-on collisions, 
in which the black holes exhibit changes in their coordi¬ 
nate location, require a suitable guess. We must empha¬ 
size that any gauge condition we introduce should not in 
principle affect the physical content in the spacetime. In 
practice, we specify gauge conditions that facilitate the 
merger and achieve long-term-stable evolutions. 

We require that the lapse function and shift vector 
resemble those of a single, boosted Kerr-Schild black hole 
in the vicinity of each singularity, namely Eqs. and 
0 respectively. Using these single, boosted black hole 
expressions, we construct the following gauge conditions: 

Q = ( 12 ) 

(13) 

Here Q represents the densitized lapse, is the analytic 
three-metric given in Eq. ®, 7 its determinant and 
^a, ^Pj are the analytic lapse and covariant shift 

of holes A and B. It is not difficult to show that the 
expressions o and d have a close limit. That is, 
at zero separation and boost velocity, they coincide with 


the densitized lapse and shift of a single, non-rotating, 
non-boosted Kerr-Schild hole of mass M = 2m. 

It is important to emphasize that the gauge conditions 
ns and m require information about the coordinate 
location of the black hole singularities and boost veloci¬ 
ties. Here, black hole coordinate location is not meant to 
be a formal statement about the precise location of the 
physical singularities, but it is merely used to represent 
the point with respect to which HU) and m are evalu¬ 
ated. Therefore, in order to use these gauge conditions, 
there remains the task of finding suitable trajectories for 
the black holes. We have found the following “Newto¬ 
nian” trajectories adequate for this purpose: 

r{t) =ro+ vot + . (14) 

Here rg and Vq are the initial singularity position and 
boost velocity, respectively. The value for qq is obtained 
iteratively using a small number of preliminary runs to 
ensure the largest gradients of the gauge fields remain 
well covered by the excision mask. During the evolution, 
the boost velocity needed to compute m and is 
obtained simply from v = df/dt. 

B. The Merger of Excision Masks 

The key point during the merger stage is, once a com¬ 
mon apparent horizon has formed and been identified, 
to replace the individual spherical excision masks with 
a single spherical one. The size of this new excision re¬ 
gion is determined by the minimum sphere containing 
the original excision masks. However, special care must 
be taken with regard to the timing of this replacement. 
One must check that the new spherical mask is fully con¬ 
tained within the apparent horizon. Typically, when the 
common horizon forms, it has an elongated or prolate 
shape. Therefore, a delay in merging the mask allows 
the common apparent horizon to become more spherical, 
thus facilitating the creation of a spherical merged mask. 
On the other hand, there is a constraint on how long this 
delay can be. Because of the extrapolations involved in 
the excision, one cannot delay the merger of the masks 
until their surfaces touch each other. Even if the common 
apparent horizon is not completely spherical, the masks 
are merged if there are not enough grid-points between 
the individual masks to carry out excision extrapolations. 
Of course, in this case one still has to continue to sat¬ 
isfy the requirement that the merged mask is contained 
within the apparent horizon. Eor the head-on collisions 
under consideration, we find that the common apparent 
horizon becomes almost spherical before the separation 
of individual masks decreases to the point of preventing 
excision extrapolations. 

Resolution also plays a factor during the mask merger. 
If the grid resolution in the vicinity of the black holes 
is too coarse, the minimum coordinate separation of the 
individual masks required by excision extrapolations will 
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be too large. In consequence, it becomes more difficult to 
find a spherical merged mask that comfortably fits within 
the common apparent horizon. For the typical resolution 
used in our simulations, h = 0.135 M, two single-hole 
masks of radius 0.405 M require a merged mask with ra¬ 
dius of 1.1475 M. That is, single-hole masks with 40.5% 
of the single-hole horizon radius mean the new excision 
radius is approximately 57% of the common horizon ra¬ 
dius. 


C. Post-Merger Ring-Down 

After the excision mask merger, the key change con¬ 
cerns the gauge conditions. Post-merger, the spacetime 
effectively represents a single ringing black hole, and the 
problem at hand is similar to that of evolving a single sta¬ 
tionary hole. In Papers 1 and 2, we have demonstrated 
how this scenario can be simulated with long-term stabil¬ 
ity using either live gauge conditions or algebraic gauge 
using a densitized lapse. The same options apply to the 
post-merger stage of a black hole collision, and we have 
found that either approach provides long-term stability. 
For the extraction of wave forms, we find a similar sit¬ 
uation: live and densitized algebraic gauge conditions 
generate waveforms that differ only within numerical er¬ 
rors. For simplicity, we use densitized algebraic gauge 
conditions (see Paper 1 and 2 for details). 

There remains the question of how the transition from 
the superposed gauge 03 and to a single black hole 
gauge takes place. As mentioned before, the gauge con¬ 
ditions and 03 have the correct close limit (i.e. 
vanishing separation and boost velocities). Starting at 
the time when the masks are merged, we match the cur¬ 
rent “black hole positions” given by 03 to decelerat¬ 
ing trajectories given by sixth-order polynomials in time. 
These new trajectory functions are such that the coordi¬ 
nate separation and boost velocities vanish within a time 
interval (typically 10 M). The order of the polynomial 
is fixed by requiring continuity of the boost velocity and 
acceleration. 


IV. RESULTS 

The spacetime of head-on collisions of non-spinning 
black holes is by construction axisymmetric. In the 
equal-mass case there is a further reflection symme¬ 
try across the plane perpendicular to the collision axis 
halfway between the holes. Although the simulations in 
this work are intrinsically 3D, we take advantage of these 
symmetries and evolve only one octant of the numerical 
domain. One of the primary goals of our work is to per¬ 
form 3D-simulations of head-on collisions with initial sep¬ 
arations never tried before. However, determining what 
constitutes a large separation has the following complica¬ 
tion. Previous work on head-on collisions has used exclu¬ 
sively Misner or Brill-Lindquist initial data, as for exam- 


TABLE I: Setup parameters for the runs. D and L are the 
coordinate and proper initial separations, vo is the boost 
velocity and ao the “acceleration” parameter to compute the 
gauge trajectories. 


Run 

D[M] 

L[M] 


ao [M-i] 

D06 

6 

5.49 

0.289 

0.040 

DIO 

10 

10.06 

0.270 

0.014 


pie the recent long-term-stable simulations of Alcubierre 
et al. [^ 13 . To our knowledge, the only evolutions of 
binary Kerr-Schild (HuMaSh) data is the grazing colli¬ 
sion of Brandt et al. |^, where the black holes had an 
initial coordinate separation of 5 M. 

We must emphasize that there exists no unambigu¬ 
ous definition of a black hole separation; also, given that 
we are using HuMaSh initial data, direct comparison of 
our initial black hole coordinate separations with those 
from Misner or Brill-Lindquist initial data is not correct. 
An operational measure of higher physical significance is 
given by the horizon-to-horizon proper separation. That 
is, the proper 3D distance horizon-to-horizon computed 
along the collision axis within the hypersurface where the 
initial data resides. One should bear in mind that this 
proper separation is slicing-dependent. Additionally, our 
proper distance separation suffers from the effects due to 
constraint violations intrinsic to the HuMaSh data. 

TableOlgives the black hole coordinate {D) and proper 
(L) separations for the two sets of runs (labeled DOG and 
DIO) performed in this work. In addition, the table gives 
the initial boost velocity (uq) as well as the “acceleration” 
parameter (oq) used to compute the trajectories required 
to construct the gauge variables [cf. Eg. |1L4|]. As a refer¬ 
ence, the “Cook-Baumgarte” ISCO corresponds 

to a proper separation between 4.8 and 4.88 M. Simi¬ 
larly, the orbit simulation by Briigmann et al. 0 started 
from a proper separation of about 9 M. 

We evolve the initial data on a set of four nested re¬ 
finement levels, with resolution and extent as given in 
Table El Except for the outermost level, grid spacing 
and extent always increase by a factor of two. We have 
found that a grid spacing significantly larger than 1.08 M 
in the outer regions gives rise to substantial reflections of 
gravitational waves at the refinement boundary. For this 
reason we have pushed the outermost level to 140.4 M 
without adding a further coarser level. 

Using the setup described above, we obtain long-term- 
stable evolutions. The three simulation stages described 
in Sec. lIIIl are illustrated in Fig.E This figure contains 
snapshots of the DIO run at times t = 0, 10.8 M and 
14 M, including the marginally trapped surfaces as cal¬ 
culated by Thornburg’s AHFinderDirect code | ^l^ . 
The upper panel shows the initial configuration with 
holes present at z = ±5 M. As mentioned above the com¬ 
putational domain only includes the hole at 5 M due to 
symmetry. During this pre-merger stage each hole is sur¬ 
rounded by its own apparent horizon (AH). In the middle 








TABLE II: Grid parameters of the refinement levels used in 
the simulations. The number of grid points N includes ghost, 
symmetry and outer boundary zones, rmax is the extent in 
the X, y and 2-directions. 


ref-level 

h [M] 

N 

rmax [M] 

1 

0.135 

67® 

8.505 

2 

0.27 

67® 

17.01 

3 

0.54 

67® 

34.02 

4 

1.08 

CO 

CO 

140.4 



FIG. 4: Snapshots of the DIO run taken at times t = 0, 10.8 
and 14 M. From top to bottom the panels show the individual 
holes with distinct apparent horizons (only one hole is visible 
due to the octant symmetry), a common apparent horizon 
forming around the individual trapped surfaces and a single 
apparent horizon surrounding a single excision region. 


panel, the holes are sufficiently close to one another that 
a common AH has formed. However, as mentioned be¬ 
fore, it is not possible yet to switch to a single excision 
region because of the elongated shape of the common ap¬ 
parent horizon. The lower panel shows the beginning of 
the post-merger stage. 



FIG. 5: The £ = 2,m = 0 Zerilli function extracted at dif¬ 
ferent radii for the D06 (upper panel) and DIO (lower panel) 
runs. The curves have been shifted in time by the difference 
in propagation time of the wave signal between the different 
radii. 


A. Waveforms 

We extract waveforms at radii 12.5, 20, 25 and SOM 
using the Zerilli-Moncrief formalism; details are given in 
App. \B In principle, it is desirable to extract gravi¬ 
tational waves at large radii so as to ensure that the 
assumptions inherent to the Zerilli-Moncrief extraction 
are satisfied. We observe two problems, however, which 
arise out of choosing excessively large radii. First, outer 
boundary effects require less time to propagate inward 
and reach the extraction radius. Second, the actual wave 
signal contained in the metric and curvature variables is 
weaker at larger radii and thus requires higher accuracy 
in the numerical evolution. In combination with the in¬ 
creasingly coarse resolution used in outer levels, this gives 
rise to a degree of distortion and noise in the extracted 
Zerilli function which is not present at smaller radii. We 
find the values used above to be a reasonable compromise 
for the post-merger stage of the head-on collision. In 
Fig .[Slwe show the resulting f = 2, m = 0 waveforms, each 
shifted in time by the amount required by the signal to 
travel from r = 12.5 M to the corresponding extraction 
radius. All curves in the figure show a significant initial 
signal up to about t « 10 M. We trace the origin of this 
signal to the fact that the two black holes are far apart in 
the early stages of the simulation and thus the assump¬ 
tion underlying the wave extraction technique, namely 
that the spacetime is that of a single black hole plus a 
perturbation, is not satisfied with sufficient accuracy dur¬ 
ing the pre-merger phase. This viewpoint is strengthened 
by the observation that the spurious signal is stronger 
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FIG. 6: The £ = 2, m — 0 component of the real part of the 
Newman-Penrose scalar ipi extracted at different radii for the 
DOG (upper panel) and DIO (lower panel) runs. The curves 
have been shifted in time by the difference in propagation 
time of the wave signal between the different radii. 


for the larger black hole separation of case DIO in the 
lower panel of the figure. We further note in this context 
that the merger of the black holes occurs at t = 5AM 
and 11.75 M respectively for cases D06 and DIO. It is at 
about these times that the shifted waveforms extracted at 
different radii start showing better agreement with each 
other in the figures, while the spurious signal present at 
earlier times converges away at larger extraction radii. 

Because the Zerilli potential becomes zero in the limit 
of infinite extraction radii, one would expect the curves 
in the individual panels of Fig.|3 to overlap exactly in 
this limit. The figure shows that this criterion is satis¬ 
fied with increasing accuracy at larger radii. Still there 
remain small effects of back-scattering of the wave signal 
at the finite radii used in the extraction. The figure also 
reveals some noise developing at late times at larger radii; 
an analysis of the timing of the noise demonstrates that 
these are boundary effects propagating inward. Below, 
we will see that these effects limit the convergence prop¬ 
erties of our code and thus the reliability of the results 
after the onset of noise in the waveforms. 

An alternative to the Zerilli-Moncrief approach to wave 
extraction is provided by the Newman-Penrose scalar 
Given the Weyl tensor Cabcd, and an appropriately cho¬ 
sen null tetrad {I, n, m, fh}, "04 = Cabcd AC ffC rep¬ 
resents a measure of the outgoing gravitational wave con¬ 
tent. We calculate Cabcd using the “electric-magnetic” 
decomposition as described in [^. For spacetimes 
perturbatively close to Kerr, the appropriate choice of 
null tetrad is the Kinnersley tetrad In practice, the 
exact Kinnersley tetrad is not available in a numerical 


simulation. Even though recent work has suggested that 
it may be possible to identify unambiguously the Kinner¬ 
sley null directions of the numerical spacetime |^, the 
resulting Weyl scalars would still differ from their Kin¬ 
nersley values in polarization and scaling. For simplicity 
we use instead a symmetric null tetrad formed from an 
orthonormalized spherical coordinate tetrad: 

r ^ ^ (u“ + e“), ^ ^ («“ - e“), 

(eg + * e;), m“ = ^(eg - 1 e%). (15) 

Here u is the future-pointing unit hypersurface normal 
and the coordinate triad vectors were orthonormalized 
via a Gram-Schmidt process in the order e^, eg, e^. A 
symmetric tetrad of this kind (though possibly orthonor¬ 
malized in a different order) is a standard choice in ap¬ 
proximately spherically symmetric spacetimes 
and is the starting point of the Lazarus procedure | 6 ^ . 
While the numerical three-metric is used for orthonor¬ 
malization, the tetrad requires no explicit knowledge of 
a background spacetime, and so the resulting expressions 
do not involve parameters such as a background mass and 
angular momentum. At large extraction radii, the lead¬ 
ing difference between these tetrad components will be 
an overall (radius-dependent) scaling factor, which is not 
relevant for our purposes. Additionally, any coordinate 
distortion should be a first- order effect, and should not 
influence the waveforms to leading order. That said, we 
cannot quantify the relative strength at early times of 
the assumptions behind Weyl extraction of this type and 
Zerilli extraction as described earlier. 

In Fig.l^we plot the £ = 2, m = 0 quadrupole moment 
of the real part of rip^. The application of the factor r 
is convenient because the resulting function will propa¬ 
gate in the limit of large radii in analogy to the Zerilli 
function. While the waveforms in Figs .[S] El exhibit good 
agreement during the merger and ring-down phase, dif¬ 
ferences are apparent during the early in-fall stage. In 
particular the Newman-Penrose waveform is less severely 
affected by the deviations from an approximately spheri¬ 
cally symmetric spacetime. While a detailed comparison 
between the two methods of gravitational wave extrac¬ 
tion is beyond the scope of this work, it will be interest¬ 
ing to analyze the observed discrepancies in more depth 
in future work. The same applies to the distortion visi¬ 
ble in the Newman-Penrose waveform for the simulation 
DIO around t = 30 M. 

In the remainder of this work we will restrict our anal¬ 
ysis to the waveforms of Fig.^obtained with the Zerilli- 
Moncrief formalism. In order to estimate the quality 
of the waveforms, next we investigate their convergence 
properties. 
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FIG. 7: Three-level convergence analysis for the I = 2, m = Q 
Zerilli function extracted at r = 7.5 M. 


B. Convergence 

The discretization of the BSSN-evolution equations in 
the Maya code is implemented with second-order accu¬ 
racy in the grid spacing h. The wave extraction, as de¬ 
scribed in Appendix \B involves numerical approxima¬ 
tions of higher-order accuracy, such as the interpolation 
of grid functions onto the sphere of extraction and the 
evaluation of the integrals over the spheres. The nu¬ 
merical error in the wave functions should therefore be 
dominated by the second-order accuracy of the evolution 
of the field equations, and we expect the wave forms to 
exhibit second-order convergence. 

In order to verify this numerically, we focus our atten¬ 
tion on the head-on collision case D06 and carry out two 
additional runs. Compared with the resolutions listed 
in Table un the two additional runs are finer by fac¬ 
tors of 1/1.5 and 1/1.5^. On the finest grid for each 
of the three runs, this leads to resolutions hi = 0.135 M, 
/i 2 = 0.09 M and = 0.06 M, respectively. Because 
the coarsest mesh in D06 has an extent of 140.2 M, im¬ 
plementing the 1/1.5 and 1/1.5^ refinements far exceeds 
the available memory of our hardware resources. For 
this reason, instead of using mesh-refinements with the 
sizes given in Table ^ we have performed the conver¬ 
gence analysis using for each run four refinement levels 
of extent r^ax = 5.4, 10.8, 21.6 and 43.2 M. The key dif¬ 
ference compared with the simulations discussed above 
is that outer boundary effects will reach the extraction 
radii at earlier times. To mitigate this effect we extract 
the Zerilli function at the smaller radius rex = 7.5 M for 
this analysis. 

In Fig.[71we show the difference in the £ = 2, m = 0 Zer¬ 
illi function obtained for the three different resolutions. 


The difference in If 20 between the medium and high reso¬ 
lution runs has been amplified by the factor 1.5^ expected 
for second-order convergence. The resulting curve agrees 
reasonably well with the difference obtained from using 
the low and medium resolution until about 35 M, indicat¬ 
ing second-order convergence. After that time, noise orig¬ 
inating from the outer boundary at 43.2 M has reached 
the extraction radius at 7.5 M and spoils the convergence 
properties. This result is confirmed by analysis at larger 
radii where second-order convergence is observed for cor¬ 
respondingly shorter periods of time. We conclude that 
our code facilitates second-order convergent wave extrac¬ 
tion until outer boundary effects have had the time to 
propagate inward toward the radius of wave extraction. 
For the simulations presented in the previous subsection, 
this result implies reliable waveforms up to a time of 
about 100 — 125 M depending on the extraction radius. 

In view of this result it appears desirable to achieve 
more satisfactory boundary treatments. Various alter¬ 
natives have been suggested to this effect in the past. 
Conceptually the most elegant approach consists in the 
incorporation of null-infinity in the numerical domain via 
the so-called Cauchy-characteristic matching. The key 
feature of this technique is the matching of a “3-1-1” for¬ 
mulation in the interior regions containing the strong- 
held sources to a characteristic formulation of the exte¬ 
rior spacetime (see e.g. |^). This facilitates a 

straightforward compactihcation of the spacetime and, 
thus, the implementation of exact boundary conditions. 
An alternative matching technique has been suggested 
in 1^. Here the weak-held exterior spacetime is de¬ 
scribed in terms of non-spherical perturbations of a 
Schwarzschild background. A diherent strategy without 
matching involves the development of alternative bound¬ 
ary conditions which by construction preserve the con¬ 
straint equations (see e. g. mil Hi). While the deriva¬ 
tion of such boundary conditions has largely been moti¬ 
vated by stability issues, it would be interesting to see 
their effect on spurious rehections as observed in our 
simulations. The implementation of either of these al¬ 
ternative boundary treatments represents a considerable 
challenge, however, and is beyond the scope of this work. 
Instead we reduce the impact of boundary ehects by plac¬ 
ing the outer boundaries far from the strong held sources. 

Finally we emphasize that in spite of the accuracy limi¬ 
tations arising from our implementation of the boundary 
conditions, we have not observed any sign of instabili¬ 
ties: all simulations presented in this work have contin¬ 
ued for many hundreds or even thousands of M gradually 
approaching a stationary spacetime containing a single 
black hole. 


C. Quasi-normal Ringing, Masses and Radiated 
Energy 

Previous studies have shown that most of the gravita¬ 
tional radiation of a head-on collision is emitted in the 
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TABLE III: Quasi-normal ringing damping a and oscillation 
frequency w measured from the extracted waveforms. The 
final column gives the percentage of the initial ADM mass 
radiated away in the form of gravitational waves as derived 
from En.ll 71 


run 

J"ex/M 

O'Madm 

wMadm 

D,,d [%] 

DOG 

12.5 

0.0917 

0.365 

0.20 

DOG 

20 

0.0921 

0.362 

0.15 

DOG 

25 

0.0923 

0.363 

0.14 

DOG 

30 

0.0947 

0.366 

0.14 

DIO 

12.5 

0.0973 

0.369 

0.37 

DIO 

20 

0.09G3 

0.363 

0.16 

DIO 

25 

0.09G0 

0.365 

0.14 

DIO 

30 

0.097G 

0.366 

0.13 



tm 


FIG. 8: ADM mass as a function of time evaluated at var¬ 
ious radii for the DOG (upper panel) and DIO (lower panel) 
simulation. 


post-merger ring-down phase, with the total radiated en¬ 
ergy typically being less than a percent of the total ini¬ 
tial mass 1^1^ 1^1^. Furthermore, the quasi-normal 
ringing of the post-merger black hole is dominated by 
the £ = 2, m = 0 mode with small contributions from 
the £ = 4, m = 0 multipole. Our simulations represent 
the first head-on collisions with wave extraction using 
Kerr-Schild type data. They are in good agreement with 
the picture described above. 

We first note that up to £ = 3 the only multipole with 
a notable wave signal is the £ = 2, m = 0 quadrupole 
shown in Fig.^l Second, we have performed least-square 
fits of the waveforms to 

f{t)=Ae~'^*cos{ujt — (j)o). (16) 

in order to calculate the ringing damping a and oscil¬ 
lation frequency ui. The results are given in Table ITTTI 


Mah[M] 



FIG. 9: The apparent horizon mass as a function of time 
calculated for the DOG (upper panel) and DIO (lower panel) 
simulation. 


The resulting damping and oscillation frequency values 
agree within about 10%. For comparison, the quasi¬ 
normal mode spectrum of the Schwarzschild black hole 
calculated by Chandrasekhar and Detweiler predicts 
aA4 = 0.0890, wAd = 0.374, where A4 denotes the mass 
associated with the black hole background spacetime. In 
Table CIII the results are reported in units of the total 
initial ADM mass Madm- The differences between these 
two masses are the radiated energy and numerical errors. 

In the course of the evolution we numerically approxi¬ 
mate the ADM mass by evaluating the ADM integral at 
finite extraction radii. The result is shown in Fig.|^for 
simulation DOG in the upper panel and DIO in the lower 
panel. The initial values agree well with the sum 'jM 
of the special relativistic masses of the two individual 
holes moving with the velocities given in Table Ql that 
is Madm = 1.045 M and 1.039 M respectively. During 
the evolution the numerical ADM mass remains within 
a few percent of the initial value but in contrast to in¬ 
tuitive expectations it increases slightly rather than de¬ 
creases after the gravitational waves have passed the ex¬ 
traction radius. Even though this increase becomes less 
pronounced at larger extraction radii we do not consider 
this calculation of the ADM mass sufficiently accurate to 
support a calculation of the total radiated gravitational 
wave energy from conservation arguments. 

In Fig. M we plot the apparent horizon masses as calcu¬ 
lated by AHFinderDirect for the two cases DOG (up¬ 
per) and DIO (lower panel). During the pre-merger stage 
there is no common apparent horizon present and we 
plot instead the sum of the horizon masses associated 
with the individual holes (dotted curves). In contrast 
the solid curve represents the horizon mass of the sin¬ 
gle post-merger hole. In the early stages after merger 
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the horizon mass shows a significant increase in time but 
quickly settles down in the final values of 1.06 M and 
1.05 M respectively. These values show good agreement 
with the values 1.045 M and 1.039 M used for the black 
hole background mass for the wave extraction. Theoret¬ 
ically one would not expect the horizon mass to decrease 
as a function of time. We attribute the decrease by about 
one percent visible after the merger in the upper panel 
of Fig. 13 to numerical inaccuracies. 

We estimate the total radiated energy from the 
Landau-Lifshitz formula (see e.g. [^) which gives the 
radiated power as 


where the constant factor arises from our definition of the 
Zerilli function in Eq. (IA6II . Ideally this function needs 
to be evaluated at null infinity. As this is not possible 
in this type of computation, we approximate the exact 
result by calculating the radiated energy at various finite 
radii. For this purpose, we have evaluated the integral 
of P from t = Q to t = 140 M — r^x', that is, we have ex¬ 
cluded late times when the extraction radius is causally 
connected to the outer boundary. The results are shown 
in the last column of Table EH and show the impact at 
small extraction radii of the spurious signal in the Zerilli 
function at early times when the two holes are still far 
apart from each other. This feature leads to an overesti¬ 
mate of the radiated energy but quickly converges away 
if we extract at larger radii. A close investigation reveals 
that this spurious contribution drops to a few percent 
of the total radiated energy at rex = SOM in the case 
DIQ and even less in the case 1106. Within these error 
bounds we thus find a radiated energy of about 0.14 and 
0.13 percent of the total ADM mass of the system. 

V. CONCLUSIONS 

We have presented the first long-term-stable binary 
black hole collisions with wave extraction using fixed 
mesh-refinement, Kerr-Schild type initial data and dy¬ 
namical excision. While the success of our dynamical 
excision technique has been demonstrated in the past in 
the case of static and dynamic single black hole space- 
times, the results in the present work demonstrate the 
applicability of our dynamical excision to binary black 
hole systems with gravitational wave generation. 

Our simulations start from initial data constructed 
via superposition of two single Kerr-Schild black holes, 
namely HuMaSh data. We have demonstrated by conver¬ 
gence analysis of the Hamiltonian and momentum con¬ 
straints that the inherent constraint violations in the Hu¬ 
MaSh data is not significant compared with the level of 
accuracy allowed by current computational resources. 

We have conceptually divided the evolutions into three 
stages: the initial in-fall, the merger and the post-merger 
ring-down of the resulting single black hole. In particular 


during the first two stages, special care must be taken to 
ensure that the excision region always be entirely con¬ 
tained inside the apparent horizon. For this purpose, we 
have calculated the apparent horizon at each time step 
and verified that no causal violation of the excision tech¬ 
nique occurs. 

We have extracted gravitational waves using the 
Zerilli-Moncrief approach applied to a black hole back¬ 
ground spacetime in Kerr-Schild coordinates. Our results 
confirm previous findings. Wave emission is dominated 
by the quasi-normal ringing of the single post-merger 
black hole. Up to ^ = 3 the only notable contribution 
is the £= 2, m = 0 quadrupole. Comparison of the 
damping and oscillation frequencies with analytic results 
shows good agreement. We do not find the accuracy in 
calculating the ADM or the apparent horizon mass suffi¬ 
cient, however, to support a calculation of the total radi¬ 
ated gravitational wave energy purely from conservation 
principles. Instead, we use the Landau-Lifshitz formula 
to calculate the radiated energy from the Zerilli function. 
We find this energy to be of the order of 0.1 % of the ADM 
mass. Our results demonstrate the need to extract grav¬ 
itational radiation at sufficiently large radii in order to 
facilitate reasonable accuracy. In particular, we observe 
a spurious wave signal in the Zerilli function caused by 
a systematic deviation during the early pre-merger stage 
of the numerical spacetime from a perturbed single black 
hole spacetime. The wave contribution due to this spu¬ 
rious feature converges away as the extraction is carried 
out at larger radii. 

We have also calculated the Newman-Penrose scalar 
i /'4 as an alternative measure of the emitted gravitational 
radiation. While the results qualitatively confirm those 
obtained using the Zerilli approach during the merger 
and ring-down phase, differences are apparent during the 
early fall-in phase. In particular we find the Newman- 
Penrose waveforms to be less severely affected by the 
deviations from an approximately spherically symmet¬ 
ric spacetime at early stages. On the other hand the 
Newman-Penrose scalar '04 shows some distortion in the 
case of larger initial separation of the holes. In future 
work we plan to investigate in more detail the causes of 
these distortions as well as compare the different wave 
extraction techniques and their discrepancies at early 
stages. 
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APPENDIX A: WAVE EXTRACTION USING 
THE ZERILLI-MONCRIEF FORMALISM 

Our extraction of gravitational waves is based on the 
gauge invariant perturbation formalism of Moncrief M- 
A more detailed description of how this facilitates the 
extraction of gravitational waves from numerical simula¬ 
tions is given in . Applications of this method can be 
found, for example, in Iz 3) [Z^ or [t^- 

The key idea is that at sufficiently large distance from 
a strong field source, the spacetime metric can be 
viewed as a spherically symmetric background plus non- 
spherical perturbations which obey the Regge-Wheeler 


equation (odd) and the Zerilli equation (even parity per¬ 
turbations) respectively. In our case the background will 
be given by that of a single, non-rotating black hole 
in Kerr-Schild coordinates. Strictly speaking, the co¬ 
ordinates are ingoing-Eddington-Finkelstein coordinates, 
namely Kerr-Schild coordinates in the case of a non¬ 
rotating black hole. The symmetry properties of an equal 
mass head-on collision imply that the odd perturbations 
vanish, so we will restrict our attention to the even parity 
case. 

Given the Kerr-Schild background metric 

g^^dx^dx'' = — -I-dr 

-P ^1 -t_ ^ 

the first-order metric perturbations, expanded in Regge- 
Wheeler tensor harmonics [l^, can be written as 

v = (A2) 

Im 

where 




hfrdgYim 

hi^d^Yi^ 


\ 

Sym 


hi^dgYi^ 

hi^d^Yi^ 



Sym 

Sym 

r2 ^ G^^dee) 

^2Qtm d^Ytra 



^ Sym 

Sym 

Sym 

sin^ 0 -b (cot 6 dg + sin"^ 

9 d(j,^')\ 

^£m ) 


Here the spherical harmonics form a com¬ 

plete orthonormal system and are eigenvectors of the 
operator — \dgg -\- cot 9dg -I- sin“^ with eigenvalue 

+ 1). A straightforward calculation shows that the 
following linear combinations of the perturbation func¬ 
tions are invariant under first-order coordinate transfor¬ 
mations 




2Mr^ 


r-2M * 


dtK^^ -b r --T77^o 

r—2M r — 2M 


* . . s 




2/il™ 
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'(■m 


- r^drG^'^ - 


2A4r^ 
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•im 


r-2M 
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(A4) 


The linearized Einstein field equations are then equiva¬ 


lent to two equations expressed in terms of these gauge 
invariant variables. The first of these equations is the 
constraint 


dr + 2 - ^Tindt 

r - 2M 


1 - — ) i2q,+iie+l)q2) 


= 0 . 

>5) 

The second equation is more conveniently written in 
terms of the Zerilli function defined in terms of the gauge 
invariant variables by 


= 


-2(r-2M) 

[£(£ -b 1) - 2]r -b 6M 


[2qi^+£i£+l)q^2n,- (A6) 


This function represents the dynamical degree of free¬ 
dom of the polar perturbations of the black hole and 
obeys the second independent combination of the Ein¬ 
stein equations 
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drr^im 4 “ q {pt^Em Em') 4 ” — 0 , 


with the potential 


e{£ + !)[£{£ + 1) - + 6[£{£ + 1) - 2pMP + 36[£(£ + 1) - 2]M^r + 72M^ 

r^{[£{£ + 1) - 2]r + 


(A7) 


(A8) 


This is the well known Zerilli equation written in Kerr- 
Schild background coordinates and represents a special 
case of the derivation of the Zerilli equation for arbi¬ 
trary spherically symmetric spacetimes in iia based on 
the gauge-invariant formalism of Gerlach and Sengupta 
The more familiar form of the Zerilli equation (see 
e. g. the recent review 


Em Em ^E^Em — 0 (-^ 9 ) 


using Schwarzschild time t* and the tortoise coordinate 
r* is directly obtained from the coordinate transforma¬ 
tion 


r-2M 

G = t — 2A1 In———, 
r - 2M 

r* = r-|-27Wln————. (AlO) 

2A4 


Here D is an arbitrary constant of dimension length and 
Vi = il- 2Mlr)Vi. 


The derivation of the Zerilli function for given values 
of £ and m from a non-linear numerical simulation is per¬ 
formed as follows. We first calculate the physical 4-metric 
from the BSSN-variables and interpolate onto a sphere of 
constant extraction radius rex- Next we transform from 
Cartesian to spherical metric components. The key step 
is to use the orthogonality relations of the tensor har¬ 
monics which enable us to separate the contributions 
of the individual multipole moments from the metric. 
A straightforward albeit lengthy calculation shows that 
the perturbation functions are related to the spherical 


4-metric components by 


Tjlm 

Hq — 

_ 

_ 


h^m _ 

= 

_ 


gttYp^dn, 

gtrYl^dn, 

grrYg^dn, 
1 


£{£ + !) 

1 

£{£+!) 


[ g,e dgYl^ + d^Yl^dn 
J Sin 0 


gre deY;^ 
1 


9r4> 
sin^ 9 


d^Yl^dn 


£{£+!)[£{£+1)-2] 


g4>4’ 


ge^ XI 


gee - I + 2 ^^dO, 

sin^^y sin0sm0 




4 - Y)_^grn 


2 

1 


gee 


sin^ 9 


Y;^dn, 


(All) 


where a superscript * denotes the complex conjugate, 


Xgm = 2{dg-cot 9)d^Yi^ (A12) 

Wim = (ds 0 - cotOde -(A13) 

\ sin fe' / 

and the integration is performed over the sphere of con¬ 
stant extraction radius. In order to also calculate the 
derivatives in Eq. lIX^ . we extract the perturbation func¬ 
tions at radii rex ~ dr, rex and rex 4- dr on every time 
slice which contains the corresponding refinement level. 
The Zerilli function then follows directly from Eqs. 

(IA 6 II . 

We have tested the wave extraction with the analytic 
solution of Teukolsky , which describes a quadrupole 
wave in linearized general relativity. The background is 
the Minkowski metric in this case which corresponds to 
the limit A4 = 0 in the relations of this section. Because 
the quadrupole £ = 2, m = 0 is the dominant multipole 
for all simulations discussed in this work, we also choose 
the Teukolsky wave with m = 0. In prescribing the free 
function F [cf. Eq. (6) in H^], we follow the suggestion 
by Eppley and superpose an ingoing and an outgoing 
















15 



FIG. 10: The numerically calculated waveforms (dashed and 
dotted curve) for the (. = 2, m = 0 Teukolsky wave are com¬ 
pared with the analytic result (solid curve). 


wave packet. This packet is of Gaussian-like shape 

F{x) = (A14) 

with X = t ^ r. The resulting analytic expression for 
the Zerilli function 4120 is rather lengthy but calculated 
straightforwardly. We evolve these data for a = 10“® 
and cr = 1 in octant symmetry on a numerical domain of 
size 20^ with two refinement levels, the fine level having 
half the extent and twice the resolution of the coarse. 
In Fig. llOl we compare the numerically computed 'I' 2 o for 
resolutions hi = 0.25 (dashed) and /i 2 = 0.125 (dotted 
curve) in the finer refinement level. We have also dou¬ 
bled the number of points on the extraction sphere from 
80 X 40 to 160 X 80 in the 6 and (f) direction. In comparison 
the solid line represents the analytic result and demon¬ 
strates the convergence of the wave-forms. A quantitative 
evaluation shows that the convergence is fourth order in 
the early and late stages of this simulation and about 
third order in between when the wave pulse reaches the 
extraction radius rex = 15. This behavior is compati¬ 
ble with the fact that we use a fourth-order integration 
scheme for the integration on the sphere and the main 
evolution code is second-order convergent. 


[1] M. Shibata, K. Taniguchi, and K. Uryu, Phys. Rev. D 

68, 084020 (2003), gr-qc/0310030. 

[2] P. Marronetti, M. D. Duez, S. L. Shapiro, and T. Baum- 
garte, Phys. Rev. Lett. 92, 141101 (2004), gr-qc/0312036. 

[3] M. Miller, P. Gressman, and W.-M. Suen, Phys. Rev. D 

69, 064026 (2004), gr-qc/0312030. 

[4] B. Briigmann, W. Tichy, and N. Jansen, Phys. Rev. Lett. 
92, 211101 (2004), gr-qc/0312112. 

[5] M. Miller, Phys. Rev. D accepted (2005), gr- 
qc/0502087. 

[6] K. R. Eppley, Ph.D. thesis, Princeton University (1975). 

[ 7 ] L. Smarr, A. Cadez, B. DeWitt, and K. Eppley, Phys. 
Rev. D 14, 2443 (1976). 

[8] L. Smarr, in Sources of Gravitational Radiation, edited 
by L. Smarr (Cambridge University Press, Cambridge, 
1979), pp. 245-274. 

[9] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 
(1995). 

[10] H. Friedrich, Class. Quantum Gravit. 13, 1451 (1996). 

[11] S. Frittelli and O. A. Reula, Phys. Rev. Lett. 76, 4667 
(1996), gr-qc/9605005. 

[12] A. Anderson and J. W. York, Jr., Phys. Rev. Lett. 82, 
4384 (1999), gr-qc/9901021. 

[13] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59, 
024007 (1999), gr-qc/9810065. 

[14] M. Alcubierre, G. Allen, B. Briigmann, E. Seidel, and 
W.-M. Suen, Phys. Rev. D 62, 124011 (2000), gr- 
qc/9908079. 

[15] L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, Phys. 
Rev. D 64, 064017 (2001), gr-qc/0105031. 

[16] G. Yoneda and H.-A. Shinkai, Phys. Rev. D 66, 124003 
(2002), gr-qc/0204002. 


[17] O. Sarbach and M. Tiglio, Phys. Rev. D 66, 064023 
(2002), gr-qc/0205086. 

[18] C. Bona, T. Ledvinka, C. Palenzuela, and M. Zacek, 
Phys. Rev. D 67, 104005 (2003), gr-qc/0302083. 

[19] M. Alcubierre and B. Briigmann, Phys. Rev. D 63, 
104006 (2001), gr-qc/0008067. 

[20] H.-J. Yo, T. W. Baumgarte, and S. L. Shapiro, Phys. 
Rev. D 66, 084026 (2002), gr-qc/0209066. 

[21] M. Alcubierre, B. Briigmann, P. Diener, M. Koppitz, 
D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 
67, 084023 (2003), gr-qc/0206072. 

[22] U. Sperhake, K. L. Smith, B. Kelly, P. Laguna, and 
D. Shoemaker, Phys. Rev. D 69, 024012 (2004), gr- 
qc/0307015. 

[23] S. Brandt, R. Correll, R. Gomez, M. Huq, P. Laguna, 
L. Lehner, P. Maronetti, R. A. Matzner, D. Neilsen, 
J. Pullin, et ah, Phys. Rev. Lett. 85, 5496 (2000), gr- 
qc/0009047. 

[24] M. Alcubierre, B. Briigmann, P. Diener, F. Herrmann, 
D. Pollney, E. Seidel, and R. Takahashi (2004), gr- 
qc/0411137. 

[25] M. Alcubierre, B. Briigmann, P. Diener, F. S. Guzman, 
1. Hawke, S. Hawley, F. Herrmann, M. Koppitz, 
D. Pollney, E. Seidel, et al. (2004), AEI-2004-115, 
gr-qc/0411149 

[26] P. Anninos, D. Hobill, E. Seidel, L. Smarr, and W.-M. 
Suen, Phys. Rev. Lett. 71, 2851 (1993), gr-qc/9309016. 

[27] P. Anninos, D. Hobill, E. Seidel, L. Smarr, and W.-M. 
Suen, Phys. Rev. D 52, 2044 (1995). 

[28] P. Anninos, R. H. Price, J. Pullin, E. Seidel, and W.-M. 
Suen, Phys. Rev. D 52, 4462 (1995), gr-qc/9505042. 

[29] J. Baker, A. Abrahams, P. Anninos, S. Brandt, R. Price, 
















16 


J. Pullin, and E. Seidel, Phys. Rev. D 55, 829 (1997), 
gr-qc/9608064. 

[30] P. Anninos and S. Brandt, Phys. Rev. Lett. 81, 508 
(1998), gr-qc/9806031. 

[31] R. H. Price and J. Pullin, Phys. Rev. Lett. 72, 3297 
(1994), gr-qc/9402039. 

[32] A. M. Abrahams and G. B. Cook, Phys. Rev. D 50, 
R2364 (1994), gr-qc/9405051. 

[33] R. J. Gleiser, C. O. Nicasio, R. H. Price, and J. Pullin, 
Phys. Rev. Lett. 77, 4483 (1996), gr-qc/9609022. 

[34] L. E. Simone, E. Poisson, and C. M. Will, Phys. Rev. D 
52, 4481 (1995), gr-qc/9506080. 

[35] J. Baker, B. Briigmann, and M. Campanelli, Class. 
Quantnm Grav. 17, L149 (2000), gr-qc/0003027. 

[36] B. Briigmann, Int. J. Mod. Phys. 8, 85 (1999), gr- 
qc/9708035. 

[37] M. Alcubierre, W. Benger, B. Biigmann, G. Lanfermann, 
and L. Nerger, Phys. Rev. Lett. 87, 271103 (2001), gr- 
qc/0012079. 

[38] D. R. Fiske, J. G. Baker, J. R. van Meter, D.-I. Choi, 
and J. M. Centrella (2005), gr-qc/0503100. 

[39] Y. Zlochower, J. G. Baker, M. Campanelli, and C. O. 
Lousto (2005), gr-qc/0505055. 

[40] D. Shoemaker, K. Smith, U. Sperhake, P. Laguna, 
E. Schnetter, and D. Fiske, Class. Quantum Grav. 20, 
3729 (2003), gr-qc/0301111. 

[41] E. Schnetter, S. H. Hawley, and 1. Hawke, Class. Quant. 
Grav. 21, 1465 (2004), gr-qc/0310042. 

[42] G. B. Cook, Living Rev. Relativity 2000-5 

(2000), [Online Article] cited on 30 Sep 2004, 

http: / / relativity.livingreviews.org/Articles/lrr-2000-5 

[43] R. A. Matzner, M. F. Huq, and D. Shoemaker, Phys. 
Rev. D 59, 024015 (1998), gr-qc/9805023. 

[44] P. Marronetti, M. Huq, P. Laguna, L. Lehner, 

R. Matzner, and D. Shoemaker, Phys. Rev. D 62, 024017 
(2000), gr-qc/0001077. 

[45] E. Bonning, P. Marronetti, D. Neilsen, and R. A. 

Matzner, Phys. Rev. D 68, 044019 (2003), gr- 

qc/0305071. 

[46] J. W. York, Jr., in Sources of Gravitational Radiation, 
edited by L. Smarr (Cambridge University Press, Cam¬ 
bridge, 1979), pp. 83-126. 

[47] C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. 
Lett. 75, 600 (1995), gr-qc/9412071. 

[48] M. Alcubierre, B. Briigmann, D. Pollney, E. Seidel, and 
R. Takanashi, Phys. Rev. D 64, 061501 (2001), AEI- 
2001-021, gr-qc/0104020 

[49] M. Alcubierre, A. Corichi, J. A. Gonzales, D. Nunez, and 
M. Salgado (2003), gr-qc/0303086. 

[50] M. Alcubierre, A. Corichi, J. Gonzalez, D. Niifiez, and 
M. Salgado, Class. Quantum Grav. 20, 3951 (2003), gr- 
qc/0303069. 

[51] C. Gundlach and J. M. Martin-Garcia, Phys. Rev. D 70, 
044032 (2004), gr-qc/0403019. 

[52] O. Sarbach, G. Calabrese, J. Pullin, and M. Tiglio, Phys. 
Rev. D 66, 064002 (2002), gr-qc/0205064. 

[53] G. Nagy, O. E. Ortiz, and O. A. Reula, Phys. Rev. D 70, 
044012 (2004). 

[54] T. W. Baumgarte, Phys. Rev. D 62, 024018 (2000), gr- 


qc/0004050. 

[55] G. B. Cook, Phys. Rev. D 50, 5025 (1994), gr- 
qc/9404043. 

[56] J. Thornburg, Phys. Rev. D 54, 4899 (1996), gr- 
qc/9508014. 

[57] J. Thornburg, Class. Quantum Grav. 21, 743 (2004), gr- 
qc/0306056. 

[58] L. Smarr, Ph.D. thesis, University of Texas at Austin 
(1975). 

[59] L. Gunnarsen, H.-A. Shinkai, and K.-I. Maeda, Class. 

Quant. Grav. 12, 133 (1995), gr-qc/9406003. 

[60] W. Kinnersley, J. Math. Phys. 10, 1195 (1969). 

[61] C. Beetle, M. Bruni, L. M. Burko, and A. Nerozzi (2004), 
gr-qc/0407012. 

[62] J. Baker, M. Campanelli, and C. Lousto, Phys. Rev. D 
65, 044001 (2002), gr-qc/0104063. 

[63] N. T. Bishop, R. Gomez, P. R. Holvorcem, R. A. 
Matzner, P. Papadopoulos, and J. Winicour, Phys. Rev. 

Lett. 76, 4303 (1996). 

[64] R. A. d’Inverno and J. A. Vickers, Phys. Rev. D 54, 4919 
(1996). 

[65] J. Winicour, Living Rev. Relativity 2001-3 (2001), 

[Article in Online Journal] cited on 23 Feb 2001, 
http://www.livingreviews.org/Articles/Volume4/2001-3winicour 

[66] A. M. Abrahams, L. Rezzolla, M. E. Rupright, A. Ander¬ 
son, P. Anninos, T. W. Baumgarte, N. T. Bishop, S. R. 

Brandt, J. C. Browne, K. Camarda, et ah, Phys. Rev. 

Lett. 80, 1812 (1998). 

[67] S. Frittelli and R. Gomez, Phys. Rev. D 69, 124020 
(2004), gr-qc/0404070. 

[68] H. Beyer and O. Sarbach, Phys. Rev. D 70, 104004 
(2004), gr-qc/0406003. 

[69] S. Chandrasekhar and S. Detweiler, Proc. R. Soc. Lond. 

A. 344, 441 (1975). 

[70] C. T. Cunningham, R. H. Price, and V. Moncrief, ApJ 
236, 674 (1980). 

[71] V. Moncrief, Ann. Phys. 88, 323 (1974). 

[72] A. M. Abrahams and C. R. Evans, Phys. Rev. D 42, 2585 
(1990). 

[73] A. Abrahams, D. Bernstein, D. Hobill, E. Seidel, and 
L. Smarr, Phys. Rev. D 45, 3544 (1992). 

[74] K. Camarda and E. Seidel, Phys. Rev. D 59, 064019 
(1999), gr-qc/9805099. 

[75] J. A. Font, T. Goodale, S. Iyer, M. Miller, L. Rezzolla, 

E. Seidel, N. Stergioulas, W.-M. Suen, and M. Tobias, 

Phys. Rev. D 65, 084024 (2002), gr-qc/0110047. 

[76] T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 
(1957). 

[77] O. Sarbach and M. Tiglio, Phys. Rev. D 64, 084016 

( 2001 ). 

[78] U. H. Gerlach and U. K. Sengupta, Phys. Rev. D 19, 

2268 (1978). 

[79] A. Nagar and L. Rezzolla (2005), gr-qc/0502064. 

[80] S. A. Teukolsky, Phys. Rev. D 26, 745 (1982). 

[81] K. Eppley, in Sources of Gravitational Radiation, edited 
by L. L. Smarr (Cambridge University Press, Cambridge, 

1979), pp. 275-291. 




